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We analyze a controversial topic about the universality class of the three-dimensio- 
nal Ising model with long-range-correlated disorder. Whereas both theoretical and 
numerical studies agree on the validity of extended Harris criterion (A. Weinrib, 
B.I. Halperin, Phys. Rev. B 27 (1983) 413) and indicate the existence of a new 
universality class, the numerical values of the critical exponents found so far differ 
essentially. To resolve this discrepancy we perform extensive Monte Carlo simula- 
tions of a 3d Ising model with non-magnetic impurities being arranged in a form 
of lines along randomly chosen axes of a lattice. The Swendsen-Wang algorithm is 
used alongside with a histogram reweighting technique and the finite-size scaling 
analysis to evaluate the values of critical exponents governing the magnetic phase 
transition. Our estimates for these exponents differ from both previous numerical 
simulations and are in favour of a non-trivial dependency of the critical exponents 
on the peculiarities of long-range correlations decay. 
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1 Introduction 



Critical properties of structurally disordered magnets remain a problem of 
great interest in condensed matter physics, as far as real magnetic crystals 
are usually non-ideal. Commonly, in the theoretical studies, as well as in the 
MC simulations, one considers point-like uncorrelated quenched non-magnetic 
impurities [1]. However, in real magnets one encounters non-idealities of struc- 
ture, which cannot be modeled by simple point-like uncorrelated defects. In- 
deed, magnetic crystals often contain defects of a more complex structure: 
linear dislocations, planar grain boundaries, three-dimensional cavities or re- 
gions of different phases, embedded in the matrix of the original crystal, as 
well as various complexes (clusters) of point-like non-magnetic impurities [2]. 
Therefore, a challenge is to offer a consistent description of the critical phe- 
nomena influenced by the presence of such complicated defects. 

Different models of structural disorder have arisen as an attempt to describe 
such defects. In this paper we concentrate on the so-called long-range-correlated 
disorder when the point-like defects are correlated and the resulting critical 
behaviour depends on the type of this correlation. Several models have been 
proposed for description of such a dependence [3-6], a subject of extensive an- 
alytical [4-14] and numerical [6, 14-17] treatment. A common outcome of the 
above studies is that although the concentration of non-magnetic impurities is 
taken to be far from the percolation threshold, in the region of weak dilution, 
the impurities make a crucial influence on an onset of ordered ferromagnetic 
phase. Given that the pure (undiluted) magnet possesses a second-order phase 
transition at certain critical temperature Tc, an influence of the weak dilution 
may range from the decrease of T^. to the changes in the universahty class 
and even to the smearing off this transition [14] . Moreover, the critical expo- 
nents governing power low scaling in the vicinity of Tc may depend on the 
parameters of impurity-impurity correlation. 

To give an example, the Harris criterion, which holds for the energy-coupled 

uncorrelated disorder [18] is modified when the disorder is long-range corre- 
lated [5,7]. In particular, when the impurity-impurity pair correlation function 
g{r) decays at large distances r according to a power law: 

g'(r) ~ 1/r", r ^ oo (1) 



the asymptotic critical exponents governing magnetic phase transition (and 
hence the universality class of the transition) do change if [5] 



pure 



< 2/a, (2) 



where u^^^^ is the correlation length critical exponent of the undiluted magnet. 
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The above condition (2) holds for a < d, d being the space (lattice) dimension. 
For a > d the usual Harris criterion [18] is recovered and condition (2) is 
substituted by i/P"" < 2/d. 

The fact, that the power of the correlation decay might be a relevant parameter 
ai a < d can be easily understood observing an asymptotics of the Fourier 
transform Qk of g{r) at small wave vector numbers k. From (1) one arrives at 
Qk ~ k"'~'^, which for a < d leads to a singular behaviour at A; — > 0. As far 
as the small k region defines the criticality, the systems with a < d are good 
candidates to manifest changes in the critical behaviour with respect to their 
undiluted counterparts. On contrary, impurity-impurity correlations at a > d 
do not produce additional singularities with respect to the uncorrelated point- 
like impurities, therefore they are referred to as the short-range correlated. In 
turn, the disorder characterized by Eq. (1) with a < d is called the long-range 
correlated. 

There are different ways to model systems with the long-range-correlated dis- 
order governed by Eq. (1). The most direct interpretation relies on the obser- 
vation that the integer a in Eq. (1) corresponds to the large r behaviour of 
the pair correlation function for the impurities in the form of points {a = d), 
hues {a = d — 1), and planes {a = d — 2) [5]. Since the last two objects extend 
in space, the impurities with a < d sometimes are called the extended ones. 
Note that the isotropic form of the pair correlation function (1) demands ran- 
dom orientation of such spatially extended objectsQ Non-integer a sometimes 
are treated in terms of a fractal dimension of impurities, see e.g. [19]. Besides 
energy-coupled disorder, the power-low correlation decay (1) is relevant for the 
thermodynamic phase transition in random field systems [20], percolation [21], 
scaling of polymer macromolecules at presence of porous medium [22,23], zero- 
temperature quantum phase transitions [24]. 

Our paper was stimulated by the observations of obvious discrepancies in the 
state-of-the-art analysis of criticality in three-dimensional Ising magnets with 
the long-range-correlated disorder governed by Eq. (1). Indeed, since for the 
pure d = 3 Ising model z/P"'''^ = 0.6304(13) [25], the long-range correlated 
disorder should change its universality class according to Eq. (2). Whereas 
both theoretical and numerical studies agree on the validity of extended Harris 
criterion (2) and bring about the new universality class [5, 12, 15, 16], the 
numerical values of the critical exponents being evaluated differ essentially. 
We list the values of the exponents found so far by different approaches in 
table 1 and refer the reader to the section 2 for a more detailed discussion of 
this issue. Here, we would like to point out that presently the results of each 
of existing analytical approaches (Refs. [5] and [12]) is confirmed by only one 



^ Anisotropic distributions of extended impurities are treated in Refs. [3,4,7-9,11, 
14]. 
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numerical simulation (Refs. [15] and [16], respectively). To resolve such a bias, 
we perform MC simulations of a = 3 Ising model with extended impurities 
and evaluate critical exponents governing ferromagnetic phase transition. As it 
will become evident from the further account, our estimates for the exponents 
differ from the results of two numerical simulations performed so far [15, 16] 
and are in favour of a non-trivial dependency of the critical exponents on the 
peculiarities of long-range correlations decay. 



Reference 


V 


7 




V 


Weinrib, Halperin, [5] 


1 


2 


1/2 





Prudnikov et al, [12] 


0.7151 


1.4449 


0.3502 


-0.0205 


Ballesteros, Parisi, [15]" 


1.012(10) 


1.980(16) 


0.528(7) 


0.043(4) 


[15]^ 


1.005(14) 


1.967(23) 


0.524(9) 


0.043(4) 


Prudnikov et al, [16]" 


0.719(22) 


1.407(24) 


0.375(45) 


0.043(93) 




0.710(10) 


1.441(15) 


0.362(20) 


-0.030(7) 



Table 1 

The critical exponents for the three-dimensional Ising model with extended impu- 
rities for a = 2. Renormahzation group calculations: first order e, (5-expansion [5]; 
two- loop massive renormahzation scheme [12]. Monte Carlo simulations: finite-size 
scaling, combination of Wolff and Swendsen-Wang algorithms, magnetic site con- 
centration p = 0.8 [15]", p = 0.65 [15]^; short-time critical dynamics with Metropolis 
algorithm, p = 0.8 [16]", finite-size scaling with Wolff algorithm, p = 0.8 [16]''. Val- 
ues of the exponents obtained from the scaling relations are shown in italic (see 
section 2 for a more detailed discussion). 

The outline of the paper is the following. In the next section we make a 
brief overview of the results of previous analysis of the 3d Ising model with 
long-range-correlated impurities paying special attention to the former MC 
simulations. Details of our MC simulations are explained in sections 3 and 4. 
There, we formulate the model and define the observables we are interested 
in. We analyze statistics of typical and rare events taking the magnetic sus- 
ceptibility as an example in section 4. The numerical values of the exponents 
are evaluated there by the finite-size scaling technique. Section 5 concludes 
our study. 



2 Overview of previous analytical and MC results 

Currently, there exist two different analytical results for the values of critical 
exponents of the 3d Ising model with long-range-correlated impurities. The 
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first one is due to Weinrib and Halperin, who formulated the model of a m- 
vector magnet with quenched impurities correlated via the power law (1) [5]. 
They used the renormalzation group technique carrying out a double expan- 
sion m. e — A — d, 5 — A — a, considering e and 5 to be of the same order of 
magnitude, and estimating values of the critical exponents to the first order 
in this expansion. Further, they conjectured the obtained first order result for 
the correlation length critical exponent 

v = 2/a (3) 



to be an exact one and to hold for any value of spin component number m 
provided that v^™^ < 2/a. By complementing Eq. (3) with the first order 
value of the pair correlation function critical exponent, = [5], the other 
critical exponents can be obtained via familiar scaling relations. For a — 2 the 
exponents are listed in table 1. 

The second theoretical estimate was obtained by Prudnikov et al. [12] in the 
field theoretical renormalization group technique by performing renormaliza- 
tion for non-zero mass at fixed space dimension d — ?>. Their two-loop cal- 
culations refined by the resummation of the series obtained bring about a 
non-trivial dependence of the critical exponents both on m and on a. We list 
their values of the exponents for m = 1 and a = 2 in the second row of table 
1. In particular, one observes that the correlation length exponent v obtained 
in Ref. [12] differs from the value predicted by (3) by the order of 25 %. As 
for the reason of such discrepancy the authors of Ref. [12] point to a higher 
order of the perturbation theory they considered together with the methods 
of series summation as well as to taking into account the graphs which are 
discarded when the £, 5-expansion is being used. 

Let us note however, the qualitative agreement between the above analytical 
results: both renormalization group treatments, Refs. [5] and [12], predict that 
the new (long-range-correlated) fixed point is stable and reachable for the 
condition considered and hence the 3d Ising model with long-range-correlated 
impurities belongs to a new universality class. These are only the numerical 
values of the exponents which call for additional verification. In such cases 
one often appeals to the numerical simulations. Indeed, two simulations were 
performed, and, strangely enough, the results again split into two groups: 
whereas the simulation of Ballesteros and Parisi [15] is strongly in favour of 
the theoretical results of Weinrib and Halperin, the simulation of Prudnikov et 
al. [16] almost exactly reproduces former theoretical results of this group [12] 
as one can see from table 1. There, we give the results for the exponents 
obtained in simulations complementing them by those that follow if one uses 
familiar scaling relations (the last are shown in italic). 
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Since our own study will rely on the simulation technique as well, we describe 
an analysis performed in Rcfs. [15,16] in more details. In Rcf. [15] the simu- 
lations were done for two different types of the long-range-correlated disorder 
referred by the authors as the Gaussian and the non-Gaussian one. In the first 
case, the point-hke defects are scattered directly on the sites of a 3d simple 
cubic lattice according to the desired decay of the pair correlation function 
(1), a = 2. In the second case, the defects form the lines directed along the 
randomly chosen axes and, as was mentioned before, their impurity-impurity 
pair correlation function should also decay at large r as g{r) ~ 1/r^). The MC 
algorithm used in Ref. [15] is a combination of a single-cluster Wolff method 
and a Swendsen-Wang algorithm. After a fixed number of a single-cluster up- 
dates one Swendsen-Wang sweep is performed. The procedure was called a 
MC step (MCS). The number of single cluster fiips was chosen such that the 
autocorrelation time te was typically 1 MCS. For thermalization of the system 
100 MCSs were performed and then various observables were measured. The 
simulation was done on the lattice of sizes L = 8, 16, 32, 64 and 128. Then 100 
MCSs have been used for measurements on A^5=20000 different samples for 
L < 64 and on A^5= 10000 samples for L = 128. The procedure seems to be 
quite safe. 

Values of the critical exponents ly and rj obtained in Ref. [15] in the case 
of Gaussian disorder (point-like long-range-correlated defects) are given in 
table 1 for the impurities concentrations p = 0.8, p = 0.65. It was stated in 
the paper, that for the non-Gaussian disorder (lines of defects) the estimates 
of the exponents are comparable with analytic calculations of Weinrib and 
Halperin [5]. 

In Ref. [16] the MC simulation of the critical behavior of the three-dimensional 
Ising model with long-range-correlated disorder at criticality was performed 
by means of the short-time dynamics and the single-cluster Wolff method. 
Randomly distributed defects had a form of lines and resembled the "non- 
Gaussian disorder" of Ref. [15]. However, in contrast to Ref. [15], a condition 
of lines mutual avoidance was implemented. According to Ref. [16], a situation 
when crossing of the lines of defects is allowed is not described by the Weinrib- 
Halperin model. In the short-time critical dynamics method, the concentration 
of spins was chosen p = 0.8 and the cubic lattices of sizes L = 16 -r- 128 were 
considered. As a MC method was used a Metropolis algorithm. Resulting 
numbers for the exponents are given in the 5th row of table 1 . 

Additional simulations in the equilibrium state were performed in Ref. [16] to 
verify the reliability of results obtained by means of the short-time critical dy- 
namics. The single-cluster Wolff algorithm was used for simulation and a finite 
size scaling for evaluation of the critical exponents. One MCS was defined as 5 
cluster fiips, lO'' MCSs were discarded for equilibration and 10^ for measure- 
ment. Disorder averaging was typically performed over 10^ samples. Again, 
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the procedure seems to be safe. The values of critical exponents obtained in 
the simulations are quoted in the 6th row of table 1. 



With the above information at hand we started our numerical analysis as 
explained in the forthcoming sections. 



3 Model and observables 

We consider a d = 3 Ising model with non-magnetic sites arranged in a form 
of randomly oriented hues (see Fig. 1). The Hamiltonian reads: 



where (ij) means the nearest neighbour summation over the sites of a s.c. 
lattice of linear size L, J > is the interaction constant, Ising spins Si = ±1, 
Q = 0; 1 is the occupation number for the i-th site and periodic boundary 
conditions are employed. Non-magnetic sites (q = 0) are quenched in a fixed 
configuration and form lines, as shown in Fig. 1. Concentration p of the mag- 
netic sites is taken to be far above the percolation threshold. Fig. 1 shows 
one possible configuration for non-magnetic impurity lines. During the sim- 
ulations, one generates a large number of such disorder reahzations. This is 
done in our study using the following procedure. The lines of impurities are 
generated one by one, each along randomly chosen axis and the final disorder 
realization is accepted with the probability 



using rejection method. Here Prcq) is a required value for the concentration of 
magnetic sites and a is the dispersion of the resulted Gaussian distribution in 
p (see Fig. 11). The average value for the impurity concentration is equal to 
1—p. In the simulations presented here we paid special attention to the width 
of this distribution and consider both broad and narrow cases, see sections 4.2 
and 4.3). 

The observables saved during each step of the MC simulations for a given 
disorder realization will be the instantaneous values for the internal energy £ 
and magnetisation A4 per spin, defined as 



(4) 




(5) 



8 — J— ^ CiCjSiS. 



(6) 
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Figure 1. A 3d Ising model with long-range-correlated disorder in the form of ran- 
domly oriented lines of non-magnetic sites. Magnetic sites are shown by discs with 
arrows, non-magnetic sites are dashed (red on-line) and do not hold arrows. 

M = ^Y.c^S,, (7) 



where number of magnetic sites is Np = pN, total number of sites is = L'^. 
Using the histogram reweighting technique [26] we compute the following ex- 
pectation values at the critical temperature region for a given disorder real- 
ization: 

(£), {\M\), {M'), {M') (8) 



as functions of an inverse temperature (3 = 1/T. In (8) the angular brackets 
(. . .) stand for the thermodynamical averaging in one sample of disorder. 

For a given disorder configuration, one can evaluate the temperature behaviour 
of the magnetic susceptibility x and magnetic cumulants U2, U4: 



In order to refer to the physical quantities, the observables are to be aver- 
aged over different disorder configurations, denoted hereafter by an overline: 
(...). Two ways of averaging can be found in the literature, which we will con- 
sider for the case of finding the maximum value for the susceptibility. At each 
disorder realization an individual curve for the susceptibility x as a function 
of the temperature T is obtained (using histogram reweighting technique). 
In the first method of averaging, depicted in Fig. 2, a the maximal values of 
each individual curves are averaged over all disorder realizations. This type 
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of averaging, referred hereafter as averaging a, will be denoted as Xmax- An 




a) Averaging a 



b) Averaging b 



Figure 2. Two ways to perform an averaging over disorder configurations. Averaging 
a: the maximal values of magnetic susceptibility are found for each disorder realiza- 
tion (black discs) and averaged over disorder configurations. The resulting averaged 
value Xmax is shown by the red square. Averaging b: the configurational averaging 
of X is performed (red dashed curve) and a maximum of the resulting curve x|max 
is found (red square). 

alternative method is shown in Fig. 2,b: at first one performs a configurational 
averaging of the susceptibility, i.e. the single averaged curve is evaluated for 
x(T). Then the maximal value for this curve, v4|max5 is found. Hereafter, this 
method will be referred as averaging b. Both averaging a and averaging b will 
be exploited in our analysis. 

To evaluate the critical exponents, a number of characteristics will be used 
with known finite size scaling: 

• averaging a: 



{\M\) 



max' 



D 



C/4 







d(3 


max 



D 



U2 



dU2 




dp 


max 



(10) 



Here, Xmax, Du^, and Djj^ are the averaged over disorder configurations 
maximal values of the magnetic susceptibility and magnetic cumulants tem- 
perature derivatives. The value (lA^Dmax calculated for the temperature 
that corresponds to the magnetic susceptibihty maximum. 
averaging b: 



X\ 



max ' 



{\M\) 



V, 



dp 



V 



U2 



dp 



'111 



V 



M 



^log(|.M|) 

dp 



a log {M^) 

dp 



:i2) 
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where 



^4 = 1- 



W2 = l- 



(13) 



■2 ' 



^\M\) 



Note, that quantities (12) are ill-defined if averaging a is considered. More- 
over, as one can see from the definitions (9), (13) there are different ways to 
define magnetic cumulants for the disordered system (see Refs. [27] for more 
discussion) . 

We performed our numerical simulations using the Swendsen-Wang cluster 
algorithm. The main reason is, that with the diluted system one often performs 
simulations at the temperature which is far from the native phase transition 
point of each particular disorder realization. As the result, the single clusters 
can be of rather small size and the Wolff one algorithm, used by us before [33] 
is found to be less effective. The following set of lattice sizes L = 6 96 (L = 
6 -f- 128 in some cases) is used and performing the finite-size-scaling analysis 
is performed. The magnetic site concentration is chosen to be p = 0.8 which 
is far from the percolation threshold and this also allows comparison with the 
previous simulations [15, 16]. Another reason for such a choice is that for the 
3d Ising model with uncorrelated impurities the correction to scaling terms 
were found to be minimal at p = 0.8 [1,28], (obviously, this does not mean 
that p = 0.8 minimizes correction-to-scaling terms at any level of impurities 
correlation). The number of samples ranges from 10^ to 10^ for all lattice 
sizes. The simulations are performed at the finite-size critical temperature 
Tc{L) which is found from the maximum of magnetic susceptibility x each 
lattice size L and the reweighting technique for the neighbouring temperatures 
is used. Note, that in Refs. [15, 16] all simulation were performed at the critical 
temperature of an infinite system. 

Similarly to our other previous papers [29, 30] we start simulation from es- 
timating the critical temperature of a finite-size system as a function of the 
lattice size L, Tc{L). For the smallest lattice size L = 6, the preliminary sim- 
ulation is performed at pTP"'''^ first (T^'^^^ is the critical temperature for the 
pure system) . The Tc{L — 6) is located then from the maximum of the sus- 
ceptibility. At the same time, we estimate various autocorrelation times to 
control the error due to time correlations. For the next larger lattice size, e.g. 
L = 8, preliminary simulation is performed at Tc{L = 6) and then Tc{L = 8) 
is located again from the maximum of the susceptibility at L = 8. The process 
is repeated for all lattice sizes. In this way we estimate Tc{L) and the energy 
autocorrelation time r^g). The results are given in table 2. One can see that 
the finite size system critical temperature obtained via averaging a differs from 
those obtained via and averaging b. The difference is not too large, the next 
simulations were performed at the critical temperature obtained by averaging 
a. 
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L 


r{£) 


Tc{L) 


Xmax 


'^Imax 





z.oUUoi 


3.796779 


3.795500 


o 



o.iUo/0 


3.812864 


3.810200 






3.821152 


3.817200 




o.ooiyi 


3.829904 


3.823000 


ib 


o.yU4Do 


3.849101 


3.842400 


on 


4.iUDoy 


3.855457 


3.849200 


ZD 


4.ZOODD 


3.863750 


3.860100 


oo 


A QCOVQ 

4.00Z (o 


3.871662 


3.869700 


A Q 


/I A AQQ 


3.876678 


3.874100 




A A^QQn 

4.4dooU 


3.881388 


3.876900 


yo 


4.zDyoo 


3.884435 


3.881200 


128 


4.21278 


3.886679 


3.883220 






3.89124 ± 0.00057 


3.89214 ± 0.00190 






3.89108 ± 0.00137 


3.89033 ± 0.00139 



Table 2 

The critical temperature Tc{L) of finite size system, obtained from maximum of 
magnetic susceptibility and autocorrelation times of absolute value of {£) in the 
Swendsen-Wang cluster sweep. 

In evaluation of the critical indices we followed the standard finite-size-scaling 
scheme as described e. g. in Ref. [26]. The estimates for the correlation length, 
magnetic susceptibility, and spontaneous magnetization critical exponents u, 
7 and /3 can be obtained from the FSS behaviour of the following quantities: 

• averaging a: 

Du,,Du,^L''\ X^-L^'\ (14) 

• averaging b: 

In Eqs. (14), (15), the correction-to-scaling terms have been omitted. Taking 
them into account the FSS behaviour of the quantity A attains the form: 



A~L'^(i + rAL-'^), 



(16) 
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where is the leading exponent, cu is the correction-to-scahng exponent, and 
is a correction-to-scahng amphtude. Hereafter, we perform the finite-size- 
scahng analysis in several ways. First, assuming the correction-to-scaling terms 
to be small we will fit the data for the observables to the leading power laws 
(14), (15) only. Besides that, we will also perform the fits with correction-to- 
scaling terms (16) as explained below. 



4 Statistics of the observables. Preliminary estimates of the critical 
exponents 

In this section, applying the above described formalism, we give preliminary 
estimates of the critical exponents. As we will see, the estimates for the same 
exponent obtained on the base of FSS of different quantities differ from each 
other. Moreover, similar differences arise when the results obtained via aver- 
aging a and averaging b are compared. We will discuss possible reasons for 
such behaviour and will analyze them further in the forthcoming section. 

The starting point of the analysis performed in this section is connected with 
the fact that at the same concentration of impurities, different samples vary 
configurationally (geometrically) as well as in their thermodynamical charac- 
teristics. In systems with uncorrelated dilution these variations are also present 
(see e.g. Ref. [31]) although, as we will see below, these are less pronounced 
than in the case of the long-range-correlated disorder. Moreover, although the 
mean concentration of magnetic sites is well-defined for an ensemble samples 
(and chosen to be p = 0.8), the concentration of magnetic sites in each sample 
vary from the mean one. Such deviations are caused by possible intersections 
of impurity lines in different samples, besides not for all lattice sizes one can 
reach given concentration of impurities by an integer number of lines. There- 
fore, to start an analysis one has to make a choice of a set of sample specifying 
allowed values of the concentration dispersion. We start an analysis by dealing 
with the set of samples that are characterized by dispersion of concentrations 
= 10-1 

4-1 Typical, average, and rare values 

The sample-to-sample variations increase with the increase of the lattice size. 
As an example, we show in Fig. 3 the values of the susceptibility for a given 

disorder realization calculated via averaging a (Xmax? Fig- 3a-3c) and averaging 
b (xlmax? Fig. 3d-3f) for three typical lattice sizes L = 10, 26, 96. Each point in 
these plots represents the data obtained for a given sample (i.e. given disorder 
realization). Solid lines are made up of the average values, where at each ji 
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the averaging is done on the interval — jj- One can see that the averages 
saturate at jl larger than two hundred samples. From this in particular one 
can conclude that considered here number of samples N^is = 1000 is quite 
sufficient for further analysis. 




200 400 600 BOO 1000 200 400 600 BOO 1000 200 400 600 800 1000 

sample # sample # sample # 



(d) (e) (f) 

Figure 3. Values of the susceptibility Xmax, calculated via averaging a ( Fig. 3a-3c) 
and xlmax' calculated via averaging b ( Fig. 3d-3f) for different disorder realizations 
(different samples). Running average over the samples is shown by the solid line. 
From left to right: L = 10, 26, 96. 



One observes that for some samples Xmax and xlmax differ significantly from 
the average values. When this is the case, we will call such an event a rare 
one. The values that are close to the maximal one will be referred to as typical 
events. We would like to stress that for averaging b for the lattice size L = 10 
most of rare events produce the values for the susceptibility below the average 
value, while for L = 96 the situation is reversed, see figure 3. At L = 26 the 
distribution of rare events is aproximately symmetric. The situation is more 
homogenous for the averaging a. 

The averaging over too small number of disorder realizations leads to typical 
(i.e. most probable) values instead of averaged ones. Indeed, as can be seen 
in Fig. 4 for L = 96, the probability distributions of Xmax and xLax possess 
long tails of rare events with larger values of susceptibility. The samples which 
correspond to the rare events give a large contribution to the average and shift 
it far from the most probable values. The value of xLax calculated with the 
complete probability distribution in shown in fig. 4 by a dashed line (blue 
on-line). Another characteristic value shown in fig. 4 by lines (solid) are the 
median values Xmed, defined as the values of x where the integrated probability 
distributions are equal to 0.5. A distance between the average and median 
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values may serve to quantify the asymmetricity of the probabihty distribution. 




Figure 4. Probability distributions of the susceptibilities Xmax; 

(Fig. 4a-4c) and 

Xlmax> (Fig- 4d-4f) for the lattice sizes L = 10,26,96 (from left to right). The fuU 
curves (blue on-line) represent the integrated distribution. Dashed vertical lines 
(blue on-line) show Xmax, and xlmax' calculated with the complete probability dis- 
tribution. Solid lines show the median. 

Figures 5 and 6 show contributions of typical and rare events to the tempera- 
ture behaviour of the observables under discussion. As an example of typical 
events, we plot in fig. 5 the magnetic susceptibility x, Binder's cummulant f/4, 
and magnetisation (|A^|) selecting the samples where susceptibilities have the 
same values as x at Pmax for the largest lattice size L = 96. Few examples of 
rare events corresponding to large values of x are shown in Fig. 6. The bold 
lines in the figures show averages over all disorder realizations (all samples). 
The thin lines (color on-line) correspond to different samples, as shown in the 
legends to the plots 5a, 6a. 




0,25711 0,25712 0,25713 0,25714 



(a) 



(b) 



Figure 5. Examples of typical events for the largest lattice size L = 96 for magnetic 
susceptibility x (a)? Binder's cummulant C/4 (b) and magnetisation (|A^|) (c). Dif- 
ferent colours correspond to different samples. The thick lines show the averages 
(calculated via averaging h) over 10^ samples. 

It follows from the above analysis that the rare events are obviously present 
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Figure 6. Examples of rare events for the largest lattice size L = 96 for magnetic 
susceptibility x (a), Binder's cummulant C/4 (b) and magnetisation (c). Dif- 

ferent colours correspond to different samples. The thick lines show the averages 
over all 10^ samples. 

in the random variables distributions and give a contribution to the average 
values of physical quantitieslZI 

These are the average values that correspond to the observables and are to 
be used for for the FSS analysis. We will perform such analysis in the next 
subsection. 

4.2 Evaluation of the critical exponents. Grand canonical disorder 

Before passing to evaluation of the critical exponents let us recall that as 
it was noted at the beginning of section 5 we perform the simulations for 
the samples with varying concentration of magnetic sites, keeping the mean 
concentration equal to p = 0.8 and limiting dispersion of concentrations by 
= 10~^ (5). Following Refs. [34] it is natural to call such situation the 
'grand canonical disorder'. Afterwards, in section 4.3 we will consider the case 
of 'canonical disorder', when we will choose the samples much more close in 
the concentration of magnetic sites (taking cr^ = 10~^). 

4-2.1 Averaging a 

In the finite-size-scaling technique [26], the critical exponents are calculated 
by fitting the data set for the observables evaluated at various lattice sizes to 
the power laws (14), (16) as was already mentioned in section 3. This pro- 
cedure will be used for the evaluation of the critical exponents u, f3, and 7 
based on the observables measured during the MC simulations. Corresponding 

^ The use of the familiar Scr law to choose the samples during the MC experiment 
(as it was done at the initial phase of our analysis, see Ref. [33]) modifies the 
distribution. Such modifying leads to different results for the exponents. We thank 
the referee for the comments on this point. 
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finite-size-scaling plots are shown in Figs. 7, 8, the numbers are given in table 
3. In Fig. 7, we give the log- log plots of the maximum values of the configura- 
tionally averaged derivatives of Binder cumulants Djj^ and Du^. The scaling of 
these quantities is governed by the correlation length critical exponent z/, see 
Eqs. (14). Fig. 8 shows the size dependence of the configurationally averaged 
magnetic susceptibility Xmax (Fig. 8a) and magnetisation (lA^Dmax (Fig- 8b). 
These values are expected to manifest a power-law scaling with exponents (3/v 
and 7/1/, see Eq. (14). 

4.5 
3.5 
2.5 
1.5 
0.5 

1.7 2.2 2.7 3.2 3.7 4.2 4.7 

InL 

Figure 7. Log-log plots for the maximum values of the configurationally averaged 
derivatives of Binder cumulants Du^ (discs) and Du^ (triangles). Solid line: data fit 
to the power law with the correction-to-scaling, Eq. (16). Here and below the range 
of the confidence interval is smaller than the symbol size. 

5.3 - 

4.3 - 

43.3 - 

^2.3 - 

1.3 - 

0.3 L 

1.5 2.0 2.5 3.0 3.5 4.0 4.5 i^j^ 1.5 2.0 2.5 3.0 3.5 4.0 4.5 

(a) (b) 

Figure 8. Log-log plots for the configurationally averaged maximum values of the 
magnetic susceptibility Xmax and magnetisation (jA^Dmax- Solid line: data fit to the 
power law with the correction-to-scaling, Eq. (16). 

We used several ways to extract the values of the exponents from the data 
given in Table 3 and plotted in Figs. 7, 8, an outcome is summarized in the 
table 4. 

First, we used the fit to the power laws (14) for all data points. Resulting 
exponents are given in the third column of table 4. Having a possibility to 
estimate v from the scaling of two different magnetic cumulants, Du^ and 
Djj^, we give in the table also an average value of V, that results from 
these estimates. However, such a straightforward estimate of the exponents 
made for the lattice sizes L = 6 -i- 96 may be not accurate enough, as far as 
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Xmax 


(l-^l)max 






6 


1.886 ±0.003 


0.4584 ± 0.0002 


8.248 ±0.007 


2.900 ± 0.003 


8 


3.312 ± 0.007 


0.3976 ± 0.0003 


12.18 ±0.02 


4.320 ± 0.007 


10 


5.068 ±0.013 


0.3533 ± 0.0003 


16.41 ± 0.04 


5.863 ±0.014 


12 


7.121 ±0.023 


0.3205 ± 0.0003 


20.50 ±0.06 


7.356 ±0.024 


16 


11.95 ±0.05 


0.2735 ± 0.0004 


28.87 ±0.13 


10.45 ±0.05 


20 


17.85 ± 0.08 


0.2413 ± 0.0004 


37.33 ±0.23 


13.65 ± 0.08 


26 


28.66 ±0.14 


0.2096 ± 0.0004 


50.10 ±0.49 


18.46 ±0.15 


35 


48.07 ±0.28 


0.1778 ±0.0005 


66.14 ±0.71 


24.58 ±0.25 


48 


86.34 ±0.55 


0.1491 ±0.0005 


92.83 ± 1.05 


34.92 ± 0.38 


64 


147.3 ± 1.0 


0.1272 ±0.0005 


125.4 ± 1.9 


48.45 ± 0.72 


96 


310.5 ±2.6 


0.1001 ±0.0004 


197.1 ±4.5 


74.71 ± 1.49 



Table 3. MC simulation data obtained for 3d Ising model with long-range correlated disorder: averaging a. 



Exponents 


Observables 


Data fit (Number of data points n 

and value of u) 


n = 11, 
oj = 


n = 5, 
u = 


n = 11, 
u = 0.8 


y 

V 




0.804(17) 
0.789(16) 


0.968(19) 
0.935(17) 


1.009(15) 
0.977(13) 


V 




0.796(12) 


0.951(13) 


0.993(10) 


15/u 


Xmax 
(l-^l)max 


1.845(11) 

0.535(5 ) 


1.825(13) 
0.560(4 ) 


1.748(16) 
0.586(5 ) 


213 /v + -ijv 


Xma^Ak^'Dmax 


2.916(15) 


2.945(15) 


2.916(19) 



Table 4 

Values of the critical exponents obtained from FSS of different observables via 
averaging a. 3rd and 5th columns: fit to all n = 11 data points, 4th column: fit to 
the 5 last data points. 

the correction-to-scahng or crossover effects may be present (for an analysis 
of crossover effects in a similar context, sec Rcf. [35]). Indeed, one can see a 
bending of the curves in Fig. 7 for L = 6^12. Taking that correction to scaling 
is especially pronounced for the small lattice sizes we decided to make a power 
law fit for five largest lattice sizes only. The results for the exponents are given 
in the fourth column of table 4. Another way to deal with the correction-to- 
scaling phenomena is to explicitly take into account the correction-to-scaling 
term in the fit via formula (16). Using theoretical estimate for the correction- 
to-scaling exponent of the 3d Ising model with long-rangc-corrclated disorder 
at the value of correlation parameter (see Eq. (1)) a = 2, u{a = 2) = 0.8 [23] 
and fitting all data points by the formula (16) we arrive at the values of the 
exponents given in the last column of table 4. An alternative way may be to 
keep the value of the exponent c<j as a fit parameter and fit all data points to 
Eq. (16) to ensure the same value of the leading exponent as those, obtained 
from the fit of data for five largest lattice sizes. As we have checked such a 
procedure leads to the values of exponent cu in the reasonable agreement with 
the theoretical estimate lo — 0.8 [23]. 

Comparing data for the exponents obtained by a power law fit (14) for five 
last data points (five largest lattice sizes) with the data obtained by a fit of all 
data points to the expression that takes into account the correction-to-scaling 
exponent, Eq. (16), we arrive to self-consistent results. Indeed, as one can see 
form the table 4, the average value of the correlation length exponent lies in 
the interval V = 0.95^0.99, the other exponents are in the range ^jv = 1.75-r 
1.83, iijv = 0.54 -=-0.58. Another check of the accuracy of the results obtained 
is the value of the combination of the exponents 2(5 /u -\- ^ /v which is to be 
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equal to three by a hyperscaling relation 2/5 + 7 = du. These value is given in 
the last row of the table. 



4-2.2 Averaging b 

Now let us carry out analysis being based on the same samples but to perform 
an averaging we will use the procedure described above as an averaging b. Re- 
sults of the analysis are summarized in tables 5, 6 and Figures 9, 10. Similarly 
as in the former subsection 4.2.1 we give in the table 5 the values of the observ- 
ables that are used in the FSS analysis. Now one can extract the correlation 
length critical exponent u from the FSS of maxima of four different quantities: 
temperature derivatives of logarithm of magnetization and of its square Vm, 
Vm^ (12) and of magnetic cumulants Vu^, Vjj^ (10). Corresponding plots are 
given in fig. 9. Again as in the former subsection we used different ways to fit 
data points to the power law dependence. On the one hand, we used simple 
power law (15), on the other hand, the correction-to-scaling was taken into 
account via formula (16). 



6.8 
5.8 
4.8 
3.8 
2.8 
1.8 
0.8 



1 1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 


1 1 1 1 


























1 1 1 1 



1.7 2.2 2.7 3.2 3.7 4.2 4.7 



luL 



Figure 9. Log-log plots for the maximum values of the configurationally averaged 
derivatives of Binder cumulants (discs), T>u^ (triangles up), Dm (triangles 
down) and D>m'^ (squares). Solid line: data fit to the power law with the correc- 
tion-to-scaling, Eq. (16). 




1.5 2.0 2.5 3.0 3.5 4.0 4.5 



Ini 




1.5 2.0 2.5 3.0 3.5 4.0 4.5 
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(a) (b) 
Figure 10. Log-log plots for the configurationally averaged maximum values of the 



magnetic susceptibility xlmax ^"^^ magnetisation 

the power law with the correction-to-scaling, Eq. (16). 



Solid line: data fit to 
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A.imax 


(lA^I) 


max 











6 


1.850 ± 





003 





4583 ± 





0006 


8.178 ± 





001 


2.850 ± 





001 


19.958 ± 





007 


33.248 ± 





013 


8 


3.205 ± 





008 





3974 ± 





0006 


11.98 ± 





01 


4.189 ± 





003 


30.269 ± 





005 


50.630 ± 





051 


10 


4.872 ± 





016 





3564 ± 





0005 


15.99 ± 





03 


5.641 ± 





012 


41.245 ± 





059 


68.993 ± 





159 


12 


6.806 ± 





026 





3245 ± 





0002 


19.73 ± 





05 


7.003 ± 





019 


52.326 ± 





102 


87.714 ± 





235 


16 


11.27 ± 





054 





2792 ± 





0002 


27.29 ± 





11 


9.807 ± 





042 


74.989 ± 





238 


125.94 ± 





47 


20 


16.65 ± 





089 





2474 ± 





0002 


35.22 ± 





19 


12.77 ± 





066 


98.707 ± 





388 


165.89 ± 





72 


26 


26.32 ± 





157 





2147 ± 





0003 


47.06 ± 





24 


17.20 ± 





10 


133.97 ± 





63 


225.22 ± 


1 


14 


35 


43.45 ± 





304 





1837 ± 





0004 


58.63 ± 





44 


22.02 ± 





18 


177.58 ± 


1 


05 


299.58 ± 


1 


75 


48 


77.51 ± 





583 





1535 ± 





0005 


84.83 ± 





08 


31.81 ± 





22 


256.23 ± 


1 


55 


431.99 ± 


2 


67 


64 


128.8 ± 


1 


073 





1314 ± 





0006 


105.8 ± 





9 


41.05 ± 





38 


335.01 ± 


2 


34 


565.50 ± 


3 


73 


96 


269.0 ± 


2 


457 





1061 ± 





0006 


178.4 ± 


1 


2 


67.75 ± 





49 


588.48 ± 


3 


01 


1012.7 ± 


4 


4 



Table 5. MC simulation data obtained for 3d Ising model with long-range correlated disorder: averaging b. 



The values of the exponent v for different fits are given in the Table 6 Having 
defined v from the scaling of four different quantities iT>UA^ T^U'z - '^m, T^w^) we 
find also an average value on i/, z7, and quote it in the table as well. Similarly, 



from the log-log dependence of x\uiax ^^"^ (l-^l) (shown in fig. 10) we 
extract values of the exponents ^/u and see Eq. (15). Again the values 
of the exponents that result from different fits are shown in Table 6. 









Data fit (Number of data points n 


Exponents 


Observables 




and value of 


a;) 




n = 11, 


n = 5, 


-1 1 

n = 11, 








a; = 


a; = 


u = 0.8 


V 






0.796(20) 


0.946(13) 


0.989(11) 


V 


Vm 




0.754(24) 


0.970(16) 


0.977(14) 


V 






0.827(20) 


1.006(43) 


1.051(18) 


V 






0.816(18) 


0.960(29) 


1.002(15) 


V 




0.798(10) 


0.971(14) 


1.005(7) 




Xlmax 




1.802(12) 


1.779(14) 


1.699(16) 


13/u 


{\M\) 


max 


0.525(4 ) 


0.540(6 ) 


0.567(5 ) 


213 /v + 7/i/ 


Xlmax' 


(1^1) 


max 


2.852(14) 


2.859(18) 


2.833(19) 



Values of the critical exponents obtained from FSS of different observables via 
averaging b. 3rd and 5th columns: fit to all n = 11 data points, 4th column: fit to 
the 5 last data points. 



As it was observed already in section 4.2.1, the strong correction-to-scaling oc- 
curs for the small lattice sizes. Similar phenomenon happens when one applies 
procedure of averaging b. An evidence may serve the bending of the curves 
in Fig. 9 for the small L. Therefore we conclude that the most reliable nu- 
merical data is obtained either from the FSS of five largest lattices or from 
all data points but with an account of the correction-to-scaling exponents. 
Prom table 6 we conclude that the exponent lays in the range V — 0.97-=- 1.01, 
7/1/ = 1.70 ^ 1.78, P/iy = 0.54 ^ 0.57. 

Tables 4 and 6 summarize results for the exponents obtained via different 
averaging procedures (a and b) and via different fitting procedures for the 
same samples. Recall that preparing the samples we have kept an average 
value of magnetic sites concentration fixed and equal to p = 0.8 with the 
dispersion = 10^^. The next step in our analysis will be to check how does 
the concentration fiuctuations infiuence the values of the (thermodynamical) 
critical exponents. 
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4-3 Evaluation of the critical exponents. Canonical disorder 



Let us consider now the situation when the dispersion of concentration p is 
much smaller as it was taken in the former subsection, that is, keeping the 
same average concentration of magnetic sites p = 0.8 let us consider much 
more narrow distribution of its values for separate samples. For the study, 
performed in this subsection we take the dispersion to be = 10^''. Doing so 
we introduce 'canonical disorder' in the spirit of Ref. [34]. As we will see below, 
this step is also can serve as an separate independent check of the values of 
the critical exponents we are interested in. 



loo- 
se - 
60- 

a 

^40- 
20- 

0- 



0.78 



/\/" 



0.79 



0.80 
P 



0.81 



0.82 



Figure 11. Distribution function p{p) of the magnetic spin concentration p for 'grand 
canonical disorder' and 'canonical disorder' in comparison, (L = 96). 

Similarly as it was described above, we perform an averaging over disorder 
in two different ways, via averaging a and averaging b. We extract the values 
of the exponents by the FSS dependence of different observables as explicitly 
given in Eqs. (14) - (16). Intermediate results and log-log plots of these ob- 
servables as functions of the lattice size are given in Appendices A and B. The 
data for the exponents are summarized in the tables 7 and 8 for averaging a 
and b, respectively. 

Comparing data of tables 4 and 6 with the data of tables 7 and 8 ones can 
see that the values of the exponents depend on the way the averaging over 
disorder was performed {averaging a and b) as well as on the way the disorder 
was prepared (grand canonical and canonical disorder). In particular, passing 
from the grand canonical disorder to the canonical disorder one observes an 
increase of 7/z/ and decrease of This difference in the exponents is more 
pronounced for the averaging a (of order of 8 %) and less pronounced for the 
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Exponents 


Observables 


Data fit (Number of data points n 

and value of a;) 


n = 12, 
w = 


n = 5, 
a; = 


n = 12, 
a; = 0.8 


V 
V 




0.841(15) 
0.817(16) 


0.887(31) 
0.904(19) 


0.989(15) 
0.975(12) 


V 




0.829(11) 


0.896(18) 


0.982(10) 


7/1/ 


Xmax 
(l-^l)max 


1.868(12) 
0.526(2 ) 


1.999(33) 
0.527(2 ) 


1.857(40) 
0.538(9 ) 


2(5 /u + 


Xmax, {\M\)^^ 


2.920(13) 


3.053(33) 


2.933(44) 



Table 7 

Values of the critical exponents obtained from FSS of different observables via 
averaging a. 3rd and 5th columns: fit to all n = 12 data points, 4th column: fit to 
the 5 last data points. 









Data fit (Number of data points n 


Exponents 


Observables 




and value of 


a;) 




n = 12, 


n = 5, 


n = 12, 








LO = 


LO = 


uj = 0.8 


V 






0.807(13) 


0.834(20) 


0.932(12) 


V 






0.764(18) 


0.866(18) 


0.956(12) 


V 


'Dua 




0.855(7 ) 


0.958(13) 


0.991(8 ) 


V 






0.841(18) 


0.921(21) 


1.031(17) 


V 




0.817(7) 


0.895(7) 


0.978(6) 




^ 1 max 




1.800(12) 


1.834(19) 


1.731(27) 


(3/v 


{\M\) 


max 


0.516(5 ) 


0.528(12) 


0.562(4 ) 


213 /v + -i/v 


xL^,{\M\) 




2.832(16) 


2.890(31) 


2.855(28) 



Values of the critical exponents obtained from FSS of different observables via 
averaging b. 3rd and 5th columns: fit to all n = 12 data points, 4th column: fit to 
the 5 last data points. 

averaging b (of order of 2-4 %). Deviation in the exponent u is of order of 5 
%. We further compare these values with the other data available and make 
some conclusions in the next section. 
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5 Conclusions and outlook 

As we explained in the introduction, our paper was inspired by existing con- 
tradictory results about the critical behaviour of the 3d Ising model with 
long-range-correlated impurities. Whereas both theoretical and MC studies 
[5, 12, 13, 15, 16] agree about the new universality class that arises in such a 
model, there exist an essential disagreement between the estimates for the 
critical exponents. For the value of the impurity correlation parameter a — 2 
they are summarized in table 1. The numerical data obtained so far split into 
two groups giving essentially different results for the exponents. To give an 
example, the predicted value for the correlation length critical exponent v de- 
viates between these groups within order of 30 % ranging from i/ = 1 [5, 15] 
to i/~ 0.71 [12,16]. 

To resolve existing contradictions, we performed MC simulations of the 3d 
Ising model with long-range-corrclatcd disorder in a form of randomly oriented 
lines of non-magnetic sites (impurity lines). Concentration of the impurities 
was taken to be 1 —p = 0.2. We used different ways to perform an averaging 
over disorder realizations, referred as averaging a and averaging b. Moreover, 
we used two different ways to prepare the samples: in one case we allowed for 
a wide distribution of concentration (grand canonical disorder) an the other 
case we made this distribution very narrow (canonical disorder). Swendsen- 
Wang MC algorithm and a finite-size-scaling analysis were applied to extract 
values of the critical exponents that govern magnetic phase transition in such 
a system. Our estimates for the critical exponents u, f3, and 7 are given in 
tables 4, 6, 7, 8. As we have discussed above, the most reliable data fits have 
been obtained by fitting the whole data for a given observable set taking into 
account the correction-to-scaling exponent or by fitting the data for the largest 
lattice sizes to the asymptotic scaling behaviour. These results are summarized 
in table 9. 

One can see from the table, that different ways to perform (and to analyze) 
the simulation lead to slightly different values of the exponents. Whereas tech- 
nical reasons for such deviations are obvious, there is no physical reason for 
them. That is all cases considered should correspond to the same asymptotic 
behaviour of physical quantities. Moreover, if in a simulation only one case 
were considered, there would be no way to see such deviation. Therefore, we 
find it reasonable to look for the averaged values of the exponents, as given in 
the last column of the table 9 

ly = 0.958 ± 0.004, -f/u = 1.809 ± 0.009, P/u = 0.551 ± 0.002. (17) 

Comparing our data quoted in table 9 with those previously obtained, table 
1, one sees that our analysis leads to the results differing from the existing 
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Exponents 


grand canonical disorder 


canonical disorder 


Averaged 
value 


averaging a 


averaging b 


average a 


average b 


n — 5, 


n = 11, 
u = 0.8 


n — 5, 
a; = 


n = 11, 
u = 0.8 


n — 5, 


n= 12, 
u = 0.8 


n — 5, 


n= 12, 
u = 0.8 


V 

7/1/ 
(3/v 


0.951(13) 
1.825(13) 
0.560(4 ) 


0.993(10) 
1.748(16) 
0.586(5 ) 


0.971(14) 
1.779(14) 
0.540(6 ) 


1.005(7 ) 
1.699(16) 
0.567(5 ) 


0.892(18) 
1.999(33) 
0.527(2 ) 


0.982(10) 
1.857(40) 
0.538(9 ) 


0.895(7 ) 
1.834(19) 
0.528(12) 


0.978(6 ) 
1.731(27) 
0.562(4 ) 


0.958(4) 
1.809(9) 
0.551(2) 


2/3/1/ + 7/1/ 


2.945(15) 


2.916(19) 


2.859(18) 


2.833(19) 


3.053(33) 


2.933(44) 


2.890(31) 


2.855(28) 


2.911(10) 



Table 9. Values of the critical exponents obtained via averaging a and averaging b for the grand canonical and for the canonical disorder. 
The last column gives averaged values. 



so far and therefore on the first sight does not make the situation more clear. 
However, we want to emphasize several items, opposing our results to previous 
theoretical and numerical estimates. 

As first discussed in Ref. [12], a difference between existing theoretical esti- 
mates for the exponents is caused by the fact that the rcnormalization group 
analysis of Weinrib and Halperin [5] was performed in the first order of the 
perturbation theory with subsequent conjecture that the first order result 
v — 2/ a is exact. Upcoming two- loop calculations at d = 3 [12] gave numeri- 
cal results for the exponents which disagree with this conjecture and question 
it. Although one may consider this result as a signal about non-trivial depen- 
dence of the exponents on the disorder correlation parameter a (as well as on 
the order parameter dimension), one certainly should not take the numbers 
obtained as a final estimates of the exponents. Indeed in the rcnormalization 
group theory of disordered systems, and, more general, of systems described 
by effective Hamiltonians of complicated symmetry, one finds many examples 
when the two-loop numerical estimates are essentially improved by higher- 
order contributions (see e.g. reviews [1] and references therein). In support to 
this suggestion let us point to the fact, that the two- loop theoretical estimates 
of Ref. [12] bring about the negative sign for the pair correlation function 
critical exponent 77 (see table 1), whereas for the 0^ theory the positiveness of 
f] follows from the Kallen-Lhemann decomposition. Therefore, our results for 
the exponents support general scenario found in Ref. [12], whilst the numerical 
discrepancy may be explained by possible changes in the theoretical estimates 
due to the high-order contributions. 

Another question concerns discrepancies between the results of numerical sim- 
ulations of the 3d Ising model with long-range correlated disorder at a = 2 
[15, 16]. In Ref. [16] the discrepancies were explained by the differences in the 
objects of simulation: in the simulations of Ballesteros et al. [15] the impu- 
rity lines were allowed to intersect opposite to the mutually avoiding impurity 
lines considered in the simulations of Prudnikov et al. [12]. One way to check 
whether the above difference in the impurity distributions causes any influence 
on the magnetic subsystem critical behaviour is to analyze an asymptotics of 
the function g{r) (1). As it follows from our analysis [32], both distributions 
lead to the power law asymptotics (1) with close values of a ~ 2. This ob- 
servation gives an evidence of the fact that both disorder distributions lead 
to the same values of the exponents. At least, if the difference exists it can 
not be seen for the system sizes considered. Moreover, if it exists it can not 
suffice to explain the observed discrepancy between simulations [15] and [16]. 
Indeed, the former simulations were performed also for the appropriately dis- 
tributed impurity sites (so-called Gaussian disorder, see section 1), when the 
very notion of the impurity lines and their intersections looses its sense. 

Let us point another phenomenon where the results obtained here may find 
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possible interpretation. For liquids, a porous medium is often better fitted to 
the extended long-range-correlated structures [36,37]. To give an example, ex- 
periments on ^He embedded in silica aerogels and xerogels report a change in 
the universality class [38] which signals about a presence of long-range corre- 
lations between defects. Simple fluids immersed in the porous media of certain 
type (determined by the density-density pair correlation function of form of 
Eq. (1)) should manifest critical behaviour in the universality class of the 3d 
Ising model with long-range-correlated disorder. Reported so far values of the 
critical exponents for the 3d Ising model with defects in a form of a porous 
medium [19, 39] were rather interpreted in terms of the random-site uncor- 
related disorder (note however a numerical agreement between the results of 
Ref. [39] and our estimates). We hope that our simulations will attract more 
attention to the interpretation of experimental and numerical studies of criti- 
cality of liquids in porous medium in terms of long-range-correlated disorder 
considered here. 

Let us return back to the values of the exponents quoted in the last column 
of table 9. From this data one gets the following values of the magnetic sus- 
ceptibility and magnetization critical exponents: 7 = 1.733(11), /3 — 0.528(3). 
Comparing these values with the results of former MC simulation of Balles- 
teros and Parisi [15] (c.f. Table 1) one can see that our result for f3 nicely 
coincide however the value for 7 differs of the order of 13 %, the difference 
in 1/ being less dramatic. Note however that our results are in much worse 
agreement with the simulations of Prudnikov et al. [16]. We do not see any 
obvious explanation for the discrepancies between MC numerical estimates of 
the critical behaviour of the 3d random Ising model with long-range-correlated 
disorder obtained so far and those obtained in our study. Eventually, if we are 
allowed to make a "not very serious" statement, we must confess that we sim- 
ply ignore the reasons of the discrepancy between previously available MC 
simulations and those presented in this work. All three studies seem equally 
rigorous, but eventually our analysis provides results for the critical exponents 
which stand between former determinations, and this may be seen as a wise 
behaviour! 

The authors acknowledge useful discussions with Christophe Chatelain and 
Wolfhard Janke. We are especially thankful to the unknown referee for illu- 
minating comments (that made our simulations even more time consuming). 
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Appendix A 



In this appendix, we give data obtained via averaging a during MC simulations 
for the canonical disorder, as explained at the beginning of section 4.3. In table 
10 we give numerical values of maxima of averaged quantities Xmax, (l-^l)max' 
Du^, and Du^, Eq. (10). Log-log plots of the dependence of these quantities 
on the lattice size L are given in Figs. 12 and 13. 
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Figure 12. Log- log plots for the maximum values of the configurationally averaged 
derivatives of Binder cumulants Du^ (discs) and Dfj^ (triangles). Solid line: data fit 
to the power law with the correction-to-scaling, Eq. (16). 
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Table 10. MC simulation data obtained for 3d Ising model with long-range correlated disorder: averaging 
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Figure 13. Log-log plots for the configurationally averaged maximum values of the 
magnetic susceptibility Xmax and magnetisation (|A^|)max- Solid line: data fit to the 
power law with the correction-to-scaling, Eq. (16). 

Appendix B 




In this appendix, we give data obtained via averaging b during MC simulations 
for the canonical disorder, as explained at the beginning of section 4.3. In table 
11 we give numerical values of maxima of averaged quantities xLax' (l-^l) ' 
T>u^, 'Du2, I^M, I^AP, Eqs. (11), (12). Log-log plots of the dependence of these 
quantities on the lattice size L are given in Figs. 14 and 15. 
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Figure 14. Log-log plots for the maximum values of the configurationally averaged 
derivatives of Binder cumulants Vjj^ (discs), Vu^ (triangles up), Vm (triangles 
down) and "Dj^.p (squares). Solid line: data fit to the power law with the correc- 
tion-to-scaling, Eq. (16). 
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Table 11. MC simulation data obtained for 3d Ising model with long-range correlated disorder: averaging b. 
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(a) (b) 
Figure 15. Log- log plots for the configurationally averaged maximum values of the 



magnetic susceptibility xlmax magnetisation (|A^|) 
the power law with the correction-to-scaling, Eq. (16). 
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